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Abstract. We present a systematic method for constructing two-dimensional magnetohydrodynamic equilibria 
with compressible flow in Cartesian geometry. This systematic method has already been developed in spherical 
geometry and applied in modelling solar and stellar winds and outflows (Vlahakis & Tsinganos, 1998) but is 
derived here in Cartesian geometry in the context of the solar atmosphere for the first time. Using the method 
we find several new classes of solutions, some of which generalise known solutions, including the Kippenhahn & 
Schluter (1957) and Hood & Anzer (1990) solar prominence models and the Tsinganos, Surlantzis & Priest (1993) 
coronal loop model with flow, and some of which are completely new. Having developed the method in full and 
summarised the several classes of solutions, we explore in a some detail one of the classes to illustrate the general 
construction method. From one of the new classes of solutions we calculate two loop-like solutions, one of which 
is the first exact two-dimensional magnetohydrodynamic equilibrium with trans- Alfvenic flow. 
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1. Introduction 

The investigation of plasma equilibria is one of the most 
important problems in magnetohydrodynamics (MHD) 
and arises in a number of fields including thermonuclear 
fusion, astrophysics and solar physics. At present there 
are many difficulties surrounding the description of fully 
three-dimensional configurations and so it is necessary to 
consider configurations with additional symmetry from a 
mathematical point of view (but see Petrie & Neukirch, 
1999). In many astrophysical situations, e.g. solar or stel- 
lar winds and outflows, axial symmetry is important while 
in solar physics translational symmetry is common in 
models of arcades, loops and prominences. 

Even with the steady-flow and symmetry assumptions, 
the prospect of solving the full MHD equations is not good 
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without further simplification. Therefore most modellers 
of MHD structures with flow adopt a numerical approach 
(e.g. Karpen et al., 2001). However we choose to use ana- 
lytical techniques as far as possible because of their greater 
computational economy and ease of use and the smaller 
role played by boundary conditions in defining a solution 
as opposed to numerical techniques, as well as the compre- 
hensiveness and clarity with which analytical methods can 
describe solution spaces. In this work we impose on our so- 
lutions spatial self-similarity i.e. the assumption of a scal- 
ing law of one of the variables on one of the coordinates. 
With this assumption the full MHD equations can be re- 
duced to a system of ordinary differential equations (ODE) 
which can be integrated by standard methods. Distinct 
from temporally self-similar solutions such as those by 
Neukirch & Cheung (2001) and the steady solutions ob- 
tained via a transformation method by Imai (1960) and 
Webb et al. (1994), many known two-dimensional equi- 
libria are spatially self-similar both in spherical geometry 
(Sauty & Tsinganos, 1994; Vlahakis & Tsinganos 1998, 
1999, Vlahakis et al. 2000) and in Cartesian geometry (de 
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Ville & Priest, 1991; Tsinganos, Surlantzis & Priest, 1993; 
Del Zanna & Hood, 1996; Del Zanna & Chiuderi, 1996). 
Two-dimensional equilibria with flow have been used in 
the context of the solar atmosphere to model coronal loops 
and arcades (de Ville & Priest, 1991; Tsinganos, Surlantzis 
& Priest, 1993) and prominences (Del Zanna & Hood, 
1996). 

A significant proportion of the energy emission from 
the solar corona is concentrated along well-defined curved 
paths called loops (Bray et al., 1991). Coronal loops are a 
feature of active regions and it is believed that they spread 
themselves out to dominate the lower corona, particularly 
in and over active regions. The loops are believed to trace 
closed lines of force of the magnetic field which penetrate 
the photosphere from below and expand to fill the whole 
of the coronal volume above an active region. Two pat- 
terns of steady plasma flow are possible in a loop: a flow 
down both legs starting at the apex, and a flow up one leg 
and down the other. Upward flows in both legs of loops 
associated with flares have also been observed but these 
are time-dependent. It has been noted that unidirectional 
loop flows are characteristic of new and complex active 
regions (Bray et al., 1991). The motion along a loop may 
continue for up to several hours so that a steady-flow de- 
scription is valid. 

In this paper we present several new classes of two- 
dimensional MHD equilibria in a uniform gravitational 
field, with flow in Cartesian geometry. The flow may be 
strictly sub-Alfvenic, super- Alfvenic or trans- Alfvenic. In 
the trans- Alfvenic case the solution can be made to cross 
the Alfvcn critical point smoothly. Discontinuities are pos- 
sible in the solution but we choose to avoid them to 
demonstrate that this is possible. We have no observed 
trans- Alfvenic loops to compare with and shocks may oc- 
cur at critical points in reality but our solution method 
allows smooth as well as shocked solutions. The solutions 
are calculated using a systematic nonlinear separation of 
variables construction method already seen in spherical 
geometry (Vlahakis & Tsinganos, 1998) and developed for 
the first time in Cartesian geometry in this paper. The 
resulting classes of solutions are found to include as spe- 
cial cases several well-known two-dimensional equilibria 
as well as many hitherto unknown solutions. To illustrate 
the general construction method we present two example 
solutions of a loop- like structure with unidirectional flow. 
We present a sub-Alfvenic case and a trans- Alfvenic case. 

The paper is organised as follows. We introduce the 
MHD equations briefly and give details of our assumptions 
in Sections 



2.1 



and 2.2 before using a simple theorem to 
construct several classes of translationally self-similar so- 
lutions in Section |2.3| , summarised in Table [l]. In Section 
H a special case of one of the classes is explored in a little 
detail by way of illustrating the general solution method 
and providing for the first time an exact loop-like solution 
with trans- Alfvenic flow as well as a new sub-Alfvenic ex- 
ample. The paper is summarised in Section ^. 



2. Construction of the model 

In this section we describe in some detail how our model 
can be systematically obtained from the closed set of gov- 
erning full MHD equations. 

2.1. Governing equations 

The dynamics of flows in solar coronal loops may be de- 
scribed to zeroth order by the well known set of the steady 
(d/dt = 0) ideal hydromagnetic equations: 



1 



p(V -V)V = — (VxB)xB 



VP-pgZ, 

V-B = 0, V-(pV) = 0, Vx(VxB)=Q : 



(1) 
(2) 



where B, V , — gZ denote the magnetic, velocity and (uni- 
form) external gravity fields while p and P the gas density 
and pressure. 

The energetics of the outflow on the other hand is gov- 
erned by the first law of thermodynamics : 



q = pV 



V 



1 p 



p. 



(3) 



where q is the net volumetric rate of some energy 
input/output (Low & Tsinganos 1986), while T — c p /c v 
with c p and c v the specific heats for an ideal gas. 

With translational Y-symmetry in Cartesian coordi- 
nates (Z, X, Y), the coordinate Y is ignorable (d/dY = 0) 
and we may introduce the poloidal magnetic flux function 
(per unit length in the Y direction) 



B = VA(Z, X) x Y + B y (Z, X)Y, 



(4) 



such that several free integrals of A exist. They are the 
ratio of the mass and magnetic fluxes on the poloidal plane 
(Z-X), V A (A), 



V - * A R 



(5) 



where the stream function ^ is a function of the magnetic 
flux function A and ^> a is its derivative, the induction 
potential $(A), V x B = V$, and the function L(A) in 
terms of which the field components in the Y-direction are 

B A 

related as V y - = L(A), and V y - —^-B y = $a(A), 
'I' \ 



(Tsinganos 1982), or, 



B 



Airp 



l-*\/^p ' 



(G) 



y -~ A l-*\/^p 

Then, the system of Eqs. (|l|) - (||) reduces to a set of two 
partial and nonlinear differential equations, i.e., the Z- 
and X- components of the momentum equation on the 
poloidal plane, 



4tt 



(£?-Vi(l- ¥\iAKn\B=V\P+^\+pgZ. (7) 
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The system of equations (Q) - (|J) - (J3J) can be solved 
to give p, P, V and B if the heating function q is known. 
Similarly, one may close this system of Eqs. ([I]) - (^) - (|J), 
if an extra functional relationship of q with the unknowns 
p, P and A exists. As an example, consider the following 
special functional relation of q with these unknowns p, P 
and A (Tsinganos, Trussoni & Sauty 1992), 



7-TP 



V-Vp, 



(8) 



where 7 < T=5/3. Then, Eq.(||) can be integrated at once 
to give the familiar polytropic relation between P and p, 



P = Q(A)p\ 



(9) 



for some function Q(A) corresponding to the enthalpy 
along a surface A — const. In this special case we can 
integrate the projection of the momentum equation along 
a stream-field line A = const on the poloidal plane to get 
the well known Bernoulli integral which subsequently can 
be combined with the component of the momentum equa- 
tion across the poloidal field lines (the transfield equation) 
to yield p and A. After finding a solution, one may go 
back to Eq.(||) and fully determine the function q(Z,X). 
It is evident that even in this special polytropic case with 
7 ^ T the heating function q (not its functional form but 
the function q(Z, X) itself) can be found only a posteriori. 
Note that for 7 = T and only then the flow is isentropic. 
Evidently, it is not possible to integrate Eq.(^) for any 
functional form of the heating function q, such as it was 
possible with the special form of the heating function 
given in Eq.(g). To proceed further then and find other 
more general solutions (effectively having a variable value 
for 7 = din P(p , A)/d\np as P is no longer integrable 
as a power of p), one may choose some other functional 
form for the heating function q and from energy conser- 
vation, Eq.(|^), derive a functional form for the pressure. 
Equivalently, one may choose a functional form for the 
pressure P and determine the volumetric rate of ther- 
mal energy a posteriori from Eq.(|^), after finding the ex- 
pressions of p, P and A which satisfy the two remaining 
components of the momentum equation. Hence, in such 
a treatment the heating sources which produce some spe- 
cific solution are not known a priori; instead, they can be 
determined only a posteriori. However, it is worth keeping 
in mind that as explained before, this situation is analo- 
gous to the more familiar constant 7 polytropic case, with 
7 7^ r. In this paper we shall follow this approach, which 
is further illustrated in the following section. 

However, even by adopting this approach, the inte- 
gration of the system of mixed elliptic/hyperbolic partial 
differential equations (|l]) - (^|) remains a nontrivial under- 
taking. Besides their nonlinearity, the difficulty is largely 
due to the fact that a physically interesting solution is 
constrained to cross some critical surfaces which are not 
known a priori but are determined simultaneously with 
the solution. 



2.2. Assumptions 

In order to construct analytically some classes of exact 
solutions for coronal loops, we shall first assume in the 
case of solutions with trans- Alfvenic flows that the mag- 
netic loop is planar, B y — 0, and there exists only a ve- 
locity component along the y-direction which is constant 
on each poloidal field line, V y = $a(A) = L(A). These as- 
sumptions are necessary for a trans- Alfvenic solution to be 
finite at an Alfven point M 2 = 1 (see Eqs. (||)). Note that 
with the assumption of a planar magnetic field the func- 
tion V y = $>a(A) = L(A) is completely decoupled from the 
equations describing the flow and is simply given in terms 
of the free integral $>a(A), or, L(A). A sub- Alfvenic or 
super-Alfvenic solution may have a non-planar magnetic 
field component as well as a non-planar velocity field com- 
ponent and be finite everywhere, provided that Eqs. (JsJ) 
are satisfied, where &a(A) and L(A) are free integrals not 
necessarily equal in solutions where Alfven points are ab- 
sent. In the present paper we will briefly discuss the possi- 
ble profiles for B y and V y for our simple examples. In the 
meantime we develop loop-like solutions by concentrating 
on the problem in the x-z plane. To proceed further we 
shall make the following two key assumptions: 

1. that the Alfven number M is solely a function of the 
dimcnsionless horizontal distance x = X/Zq, i.e., 



M = 



Bl 47T/9 



(10) 



and 



2. that the poloidal velocity and magnetic fields have an 
exponential dependence on z = Z/Zq, 



A = ZqBqA (a) , a — G(x) exp {—z) , 



(11) 



for some function G(x), where Zq and Bq are constants. 
With this formulation the magnetic field has the form 



B = B aA (a) X + F{x)Z 



where 



F(x) 



1 dG{x) 
G(x) dx 



dZ 
dX 



(12) 



(13) 



is the slope of the field lines. 

We expect that assumptions (1) - (2) are physically 
reasonable for describing the properties of loop flows, at 
least close to the solar surface. 

Instead of using the free functions of a, (A, ^a), we found 
it more convenient to work instead with the two dimen- 
sionless functions of a, (gi , $2)5 



9i («) 



ctA 2 da , 



gZ [* 2 A 
52(a) = ~Bj ~ da - 



(14) 



(15) 
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Also, we shall indicate by II the dimensionless total The above equation can be put to the form 
pressure, 

B 2 



n 



8tt 



P 



8tt 



(16) 



Using the fact that A is constant along field lines we Y\ = 1 , Y 2 = g% , Y3 
can collect the inertial term and the magnetic tension 



' V ) B _ p ry . V )y = J_ t B . v)(1 _ M 2 )B 

47T 



y 1 (a)X 1 (x)+r 2 (a)X 2 (a;) + ... + y5(«)^5(x) =0, (31) 
where, 

-ag 1 , F 4 = 92 , Y 5 = -ag 2 , (32) 

and 

X% = f ,X 2 = f 2 , X 3 = f 4 ,X A = f 3 , X 5 = f 5 . (33) 



-in 
4irZ 



<*9i 



-M X 



(1 - M 2 )F - FM 2 



We have therefore reduced the MHD equations to a 
' ,(17) single differential equation whose terms all have separated 
coordinate dependence. Note that the physical quantities 



where M 2 ' = dM 2 (x)/dx and F = dF(x)/dx. The pre- can be recovered via 
vious momentum-balance equation can be written as the 
following system for II 

1 dU(x, z 
2 



dx 
ldU(x,z) 
2 dz 



otg x 



a 9l M 2 , 

(1 - M 2 )F 



(18) 



FM 



a 92jj2- ( 19 ) 



For any differentiable function T>(x, z) the following 
relations hold 



dV(x,z) dT>(x,a) 



dx 
dV(x, z) 
dz 



dx 

dV(x, a 
da 



Fa 



dV(x, a) 
da 



(20) 
(21) 



Hence, the following differential equations are obtained 
by using the coordinates x and a 



1 dU{x,a) 

2 dx 
ldU(x,a) 
2 da 



ag'xia)/^) + ag' 2 (a)f 5 (x), 
g[{a)h(x)+g' 2 (a)f 3 (x). 



in terms of the functions of x 
1 



h{x 



(i + ^ 2 ), 

2 ' F(l-M 2 ), 



2 

FM 
1 

M 2 

-(1 + F 2 )M 

F 
~M2 



2' 



(1 - M 2 )FF 



We are now in a position to integrate Eq. ( |23| ) to obtain 
an expression for the total pressure 



(29) 



U(x,a) = 2[f 2 (x) gi (a) + f 3 (x)g 2 (a) + f (x)], 



where fo(x) is arbitrary, and to substitute this into ( |22| 
to obtain 



.4 



ZqBq / \ —da, 
J V a 

Bn 



"5 2 , 



P 

P(x, a) 



Bj ag 2 
AngZo M 2 ' 

Jo - fiag 1 + f 2 gi + hg2 



El 

47T 



B = B Q \/ag' 1 X + F(x)Z 
V 



B V Y, 



'gZ QX \^M A 

92 



X + F(x)Z 



VyY. 



(34) 
(35) 

(36) 

(37) 
(38) 

(39) 



Using our assumptions Eqs. ( [Toj , |Tl| ) we have simplified 
the full MHD equations sufficiently to undertake the task 
of solving them systematically. The system of Eqs. (Q, 
(22) has been reduced to the single Eq. (^Tj). We undertake this 
task in the next section. 

(23) 

2.3. Construction of models 



Eq.(|30|) can be solved with the help of the following the- 
orem (Vlahakis 1998). 
(24) Theorem: Let F n (a), Yi{a) and Xi(x), i = 1,2, ... ,n be 
functions of a and x such that the relation 

(25) 

(26) F "( a ) = Y 1 (a)X 1 (x) + Y 2 (a)X 2 (x) + . . . + Y n (a)X n (x)(40) 

(27) holds. Then there exist constants c±, c 2 , . . . , c n such that 

( 28 ) F n (a)=c 1 Y 1 (a) + C2Y 2 (a) + ... + c n Y n {a). (41) 
Note that Eq.(|30|) is a relation of the form 

Y n (a)X n {x) + ... + Y 1 (a)X 1 (x) = , (42) 



with n=5. Either (i) X n [x) = for every x in which case 
(indicated by the digit "0" ) we have 



fd + 9i f 2 - a 9if<i + 92 j % - ag' 2 f 5 = . 



(30) Y n _ 1 (a)X n _ 1 (x) + ... + Y 1 (a)X 1 (x)=0, 



(43) 
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or (ii) X n (x) ^ in which case (indicated by the digit 
"1") we have 

,, , » X\{x) X n ^i(x) 
X n (x) X n (x) 

Therefore according to the theorem there exist con- 
stants , i = 1, . . . ,n — 1 such that Y n (a) = 

n-l 

li^Y^a) giving a condition between the functions of 

»=i 

a. Substituting in Eq.([l4|) we obtain 



J2[Xi(x) + Vi n) Xn(x)} Yi(a) = 0. 



(45) 



i=i 



In both cases (i) and (ii) we obtain a sum with n — 1 
terms and we can continue to apply this algorithm un- 
til we have only one term. Since for each of the n terms 



in Eq.(42) there are the above two possibilities we may 
obtain 2™ cases each corresponding to a number xx . . . xx 
where each x is either or 1. The number of l's is the num- 
ber of conditions between functions of a while the number 
of 0's is the number of conditions between functions of x. 

From our relation (|30| ) we can obtain 2 5 solutions, each 
corresponding to a number xxxxx. However many of them 
are not of interest to us at present. In the first applica- 
tion of the algorithm X§ = f$ = —F/M 2 which must be 
non-zero if we are to have curved field lines, and so any 
solution of interest will have the first digit 1. Meanwhile 
Y\ = 1 so that the last digit must be 0. Furthermore 
any case \xx\Q has g\ (a) constant giving B = V = 0. 
This leaves us with the four cases lxxOO. There are three 
unknown functions of x: G(x), M(x) and fo(x), and two 
unknown functions of a: g%(a) and 52(a)- The cases 10100 
and 11000 have three equations in x and two in a and 
are well-determined. The case 10000 has four mutually- 
incompatible x-dependent equations and we must dismiss 
this case, but in the case 11100, although the three equa- 
tions in a over-determine the a-dependent problem, useful 
solutions are available. We describe the a-dependent sets 
for the families 10100, 11000 and 11100 below. 

For the case 10100: 

— The first digit is 1 , so we get Y5 as a linear combination 
of Y 1 ,Y 2 ,Y 3 , Y A . 

— The second digit is 0, so the corresponding equation is 
one between functions of x. 

— The third digit is 1, so we get Y3 as a combination of 

— The two last digits give equations between functions 
of x. 

We may neglect the Y3 term in the linear combination 
of Y\ , I2 i ^3 1 Y4 which give Y5 , since Y3 itself is a linear 
combination of Yy , Yi . 



1 Note that when we reach this step, Eq. (Q0[) is of the form 
Ti.s(x)Ys+'H2[x)Y2+'Hi(x)Yi = 0, where Hi are combinations 
oiXi. 



Using Eq.(32) we get the a-dependent set for the case 
10100: 

ag' 2 = ci + c 2 gi + c 4 g 2 and ag[ = c 5 + c 6 gi . 

In a similar way we may find the a-dependent sets for 
the cases 11000 and 11100. 

We summarize the systems of ordinary differential 
equations - with respect to a - associated with the three 
families below. 



Family 



10100 



11000 



11100 



a— dependent set 

ag' 2 = a + c 2 gi + c 4 g 2 
ag'i =c 5 + cegx 

ag' 2 — a + c 2 gi + c 3 ag[ 
.92 = c 5 + c e gi + c 7 ag[ 

ag' 2 =ci+ c 2 gi 
.92 = c 5 + c 6 gi 

ag'l = C 8 + CgCf! 



(46) 



The goal of the algorithm is to find all sets of the free 
integrals in order for the a , x variables to be separable. 
If we know the free integrals as functions of a, then it is 
trivial to find the ODE's for the functions of x; we need to 
simply substitute g\ , 172 in Eq. (^0|) and equate with zero 
the coefficients of the independent functions of a in the 
resulting sum. But before performing this step it is needed 
to clarify which are the possible forms of the free integrals; 
the solutions of Eq.(^6|). 

For the family 10100: For ce ^ Owe have g\ — — c$ /cq+ 

Ca ce , while for ce = ,gi — C5 In a + C. In both cases we 

fi A ca [ c i + c 2.9i j 
can mid q 2 = a 4 / — da. 

J a C4+1 

Exploring this last equation we conclude that the case 
10100 has the following subcases (in the final expressions 
we may ignore additive constants since the integrals de- 
pend only on the derivatives of g\ , g 2 while in the expres- 
sion for the pressure, Eq.(j3?]), the unknown /o may include 
the parts associated with these constants): 



1. If c 4 7^ 0,c 6 7^ 0,c 4 7^ c 6 , 

2. If c 4 7^ , c 6 7^ , c 4 = c 6 , 

3. If c 4 = , c 6 7^ , 

4. If c 4 7^ , c 6 = , 



5. If 



c 4 = , c 6 = , 



,9i 


= -C5/C6 


+ Ca C6 


.92 


= Ca Ci + 


C2 Q ,C fi 
C 6 — C4 


.91 


= -cs/c 6 


+ Ca Ce 


.92 


= Ca Ci + 


c 2 a Ci In a 


,9i 


= -C5/C6 


+ Ca Ce 


.92 


= ci In a - 


f- ^a C6 


,9i 


— C5 In a - 




.92 


= Ca c ± - 


^lna 


,9i 


— C5 In a - 




.92 


= Ci In a - 


f ^(lna) 2 



Next we explore the 11000 and 11100 cases. Each one 
has many subcases (as case 10100 has the above five). We 
may simplify the final expressions noting that 

— Some subcase from one case may be the same with 
some other from another case. 
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— any additive constant in g\ , g 2 can be neglected; it is 
equivalent with a rename of the unknown function fo 
in Eq.©. 

— Since a = G(x) exp (—Z/Zq) we may replace a power 
of a with a; this is equivalent with a rescaling of Zq. So 
we are free to choose the exponent in a possible term 
a x . Of course if there are more than one terms of this 
kind we can choose only one exponent. 

Our final results for the free integrals and the corre- 
sponding ODE with respect to x are summarised in Table 
1. The expressions for the pressure in each case are also 
given. R 



2.4. The subcases of some known solutions 

As mentioned in Section |l| the systematic solution method 
described in the previous subsections is capable of recov- 
ering several well-known solutions as special cases of so- 
lutions of Table [l] as well as new families first presented 
in this paper. In this section we give some details of how 
this work generalises certain specific solutions well-known 
in solar physics. 

The first family of Table 1 is simplified if we choose 
C 2 = D 2 = fo = 0. Then there is only one rela- 
tion between the unknown functions G(x) , M(x) and we 
may assume a polytropic relationship between pressure 
and density. Thus, by imposing as a second relationship 
Pi = D M~ 2 i, and since l/M 2 = Ai:p/^ 2 A , ^\ oc a 2 , it 

follows that P = ^D a 2 M- 2 "i = Q{A)pi . The isother- 
mal subcase (7 = 1) has been analyzed in Tsinganos, 
Surlantzis and Priest (1993). Besides that special subcase, 
we recover the general non-isothermal polytropic case as 
a subcase of the first family in Table 1. 

The "given temperature profile" solutions of Hood & 
Anzer (1990) belong to the first family of Table 1 also, 
with C2 = D2 = fo = 0. The second relation between 
G(x) , M(x) is the given temperature profile (the prod- 
uct P1M 2 , which is proportional to the temperature, is a 
given function of x). This is also the case with Del Zanna 
& Hood (1996), whose generalisation of Hood & Anzer's 
model to include compressible flow is similar in method to 
Tsinganos, Surlantzis and Priest (1993) and is an identi- 
cal special case of our scheme with this different imposed 
closure condition. 

The classical Kippenhahn & Schliiter (1957), 
"prominence" -like solution belongs to the second family 
of Table 1, assuming further C\ = D\ = P2 = 0. If this 
is the case, A oc In a , $>a = const, while if this constant 



2 There is another way to find all possible cases of 
Table 1. Our goal is to find three ODE's (for the functions 
fo,M,G) from Eq.@. So, all possible cases of the Table 
are answers to the question: Find three independent vectors 
[111(a) , 112(a) . un(a)] in the space of functions of a, such that 
all the Yi of Eq.(|32|) belong to the resulting subspace i.e. they 
are linear combinations of 111,112,1(3. We may choose ui (a) = 1 
(since unity is one of the Yi's) and 112(a) = gi since gi cannot 
be a constant and thus it is independent from ui). 



2 




Fig. 1. Sample magnetic and flow field lines for the sub- 
Alfvenic example. All field lines are of a loop-like shape 
and, because of the self-similarity of the solution, the lines 
differ from each other only by an additive constant. This 
solution has physical limits on the A-axis defined by a 
limit imposed on the plasma pressure and density (plotted 
in Fig. ||) so that they do not become unphysically large. 
The horizontal and vertical coordinates are given in units 
of Zq, x — X/Zq and z = Z/Zq. 

is zero we find static solutions. Again there is only one 
relation between the unknown functions G(x) , M(x) and 
we can add a polytropic relationship between pressure 
and density as the second one. 

3. A new class of solutions 

It seems that a particularly interesting new case is the 
first family in Table 1. The non-polytropic subcase, | C 2 1 + 
\D 2 \ ^ 0, has two "scales": Z and 2Z n /X. In this section 
we examine this case, with C2 = 0. 
The corresponding integrals are 



gi = C\a , 52 = Dice + D 2 a 



(47) 



.4 



gZa 



{2D ia 2 + XD 2 a x }48) 



Note that we may choose Ci such that Bo is a component 
of the magnetic field at a reference point. The form of the 
pressure is 



P 
where 



Pi a 2 +P 2 a* 



Pi = Ci 

El 

M 2 



FM 2 ' -F'(l- M 2 ) - F 2 - 1 
Pa = fo ■ 



Pi = 



Di_ 

M 2 ' 
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Table 1. Translationally self-similar exact MHD solutions 



a — integrals 



ODE 



pressure 



4 77 

Po 

Pi 

Pz 



Po + Pi a 2 + P 2 a A 



gi — C\ot + C 2 a 
g 2 = P>\o? + D 2 a* 



Cife - 2Ci/ 4 + -D1/3 - 2Z?i/ 5 = 
C2/2 - AC2/4 + £2/3 - A-D2/5 = 
fo = 



47T 

Po 

Pi 
P2 



= /o 

= C 1 f 2 + D 1 f 3 
= C 2 f 2 + D2/3 



2Ci/i 
AC2/1 



gi = C\a + C 2 In a 
<72 = Pi a 2 + P2 In a 



Cxfc - 2Ci/ 4 + Pi /a - 2Pi/ 5 = 

C 2 / 2 +D 2 /^=0 

/ - C2/4 - P2/5 = 



P = P + Pia' ! + P2lna 

= /o - C2/1 

= Ci/ 2 +Pi/3-2Ci/i 

— C2/2 + ^2/3 



47T- 

Po 
Pi 
P2 



31 = Ci In a + C 2 (In a) 2 
g 2 — Pi In a + P2 (In a) 2 



Ci/J-2Ca/* + Di/3- 

C 2 / 2 +P2/3=0 

/ - C1/4 - P1/5 = 



2P2/5 = 



47T 

Po 

Pi 
p 2 



P = P + Pi In a + P 2 (In a) 

= /o - C1/1 
= Ci/ 2 + Pi/ 3 -2(7 2 /i 
= C2/2 + P2/3 



= Po + Pia 2 + P2Q 2 In a 
/o 

Ci/2+Pi/3-(2Ci + C 2 )/i 
C2/2 + P2/3 - 2C2/1 



g\ — C\a + C 2 a In a 
<?2 = Pi a 2 + P 2 a: 2 In a 



Ci/ 2 - (2d + C 2 ) / 4 + Dift - (2Pi + Da) f 5 

C 2 f 2 - 2C 2 U + Duf's ~ 2P2 fs=0 

fo=0 



= 



2 - 




Fig. 2. Sample magnetic and flow field lines for the trans- 
Alfvenic example. All field lines are of a loop-like shape 
and, because of the self-similarity of the solution, the lines 
differ from each other only by an additive constant. This 
solution has physical limits on the X-axis defined by the 
condition that the derivative F'{x) of the slope of a line 
must be negative. At the foot points of the largest loop 
this derivative is zero and the line is exactly straight but 
reverse curvature is avoided. The horizontal and vertical 
coordinates are given in units of Zq, x — X/Zq and z = 
Z/Z . 



Using the above definitions for the "pressure components" 
together with the ODE's from Table 1, we conclude that 
for the nontrivial case A ^ , A ^ 2 , D2 ^ , we have the 
following system of equations for the six unknown func- 
tions of x (including the slope of the field lines F): 



Po 

Pi 

dlnG 

dx 
d\nM 2 
dx 

F' 



const , (49) 

IF* (50) 

F, (51) 
AP , or, M 2 = D 3 G X , with D 3 = const , (52) 



(AM 2 - 1 P 2 - 1 



Di Pi 



C1M 2 Ci 



1 — M 2 



(53) 



P[ = 



2D X F 
M 2 



- IdFM 2 [F' + X(F 2 + 1)] . (54) 



Note that practically the solution for the shape of the 
loop is a three-parametric family depending on the con- 
stants D1/C1, A and D 3 . These ODE's can be solved 



by standard numerical methods. Eqs. ( pl| , p3 54) form a 
closed system of ODE's in the three variables G, F and 
Pi. On solving these, M 2 and P2 can be recovered from 
the remaining equations. We remark that in the case 
Dx/Ci 7^ 0, F' and P[ become large for small M and the 
solution becomes unphysical. In order to obtain a loop- like 
solution it is necessary to keep F' negative. If the solution 
is super-Alfvenic near the apex of the loop where F rs 
then F' can only be kept negative by having D\jC\ 7^ 0. 
A strictly sub-Alfvenic solution can, however, have a loop- 
like appearance with D\jC\ = as will be seen, although 
we cannot save the pressure and density from becoming 
infinite as x — > ±00. From Eq.(|52|) with A > the Alfven 
Mach number graph resembles a typical loop-like field line, 
an inverted U symmetric about x = 0, but always strictly 
positive. From Eq.([l0|) p oc 1/M 2 along a field line and 
therefore p tends to 00 as x — > ±00, as does P2 from 
Eq.(|o|). In the sub-Alfvenic case it is possible to control 
P2's contribution to the solution by sending D2 — * but 
there is little advantage to be gained from doing this. Pi 




1 1 1 1 1 1 c — — 1 1 :: =- J 3 

2 4 2 4 

s/Z s/Z 

Fig. 3. Graphs of plasma quantities for the sub-Alfvenic example: plasma and magnetic pressures (top left), plasma 
density (top right), plasma temperature (bottom left) and plasma velocity and Alfven speed (bottom right) along the 
largest loop shown in Fig. [|. All quantities are graphed against the arc length s of this loop, taking as s = the left 
foot point. The U-shape profiles for the density and pressure reflects the fact that the density and pressure decrease 
with height, while the temperature and velocity increase with height along the loop. Note that the velocity is clearly 
sub-Alfvenic and that the magnetic pressure clearly exceeds the plasma pressure. 



also becomes infinite because of the dependence of P[ on 
F, which becomes infinite if F' is kept negative. 



We would like the Alfven speed Va = 
increase with height and the plasma j3 = 87rP/|B| 2 to 
decrease with height, as is the case in the solar atmo- 
sphere, and this can be guaranteed by setting A > 2. We 
remark that an alternative approach is to set < A < 2 
and D 2 < 0. We refrain from doing this since using a 
negative pressure "component" introduces complications 



in ensuring that the pressure is everywhere non-negative 
and decreasing with height. ^] 

We generate loop-like solutions as follows. We begin 
by calculating the right half of the loop, beginning from 
the loop apex at x = with suitable initial conditions 
F = 0, G > and Pi > 0. For simplicity we assume 
throughout that C\ = 1.0. The symmetry properties of 
Eqs. ([l9] - |54|) ensure that on integrating from x = 
in the negative direction the other half of a symmetric 
loop-like solution is obtained. In the sub-Alfvenic case the 



3 the above requirements are equivalent to A > and (A — 
2)I?2 > 0. So, there are two possible cases: A > 2 , D2 > or 
< A < 2,D 2 < 0. 




Fig. 4. Graphs of plasma quantities for the trans- Alfvenic example: plasma and magnetic pressures (top left), plasma 
density (top right), plasma temperature (bottom left) and plasma velocity and Alfven speed (bottom right) along the 
largest loop shown in Fig. ||. All quantities are graphed against the arc length s of this loop, taking as s — the left 
foot point. The U-shape profiles for the density and pressure reflects the fact that the density and pressure decrease 
with height, while the velocity increases with height along the loop. Note that the magnetic pressure is greater than 
the plasma pressure, and the velocity is clearly trans- Alfvenic, being super- Alfvenic near the apex and sub-Alfvenic 
near the foot points, and that the critical points are crossed smoothly. In contrast to the sub-Alfvenic example of Fig. 
U the temperature decreases with height along the super- Alfvenic part of the loop. 



equations have no critical points and can be integrated 
without difficulty. Figure [l] shows some example field lines 
of our sub-Alfvenic example and the Alfven Mach num- 
ber profile. The parameter values used in this example are 
Di/d = 0.0, D 2 = 1.0, D 3 = 1.0 and A = 3.0. Figure | 
shows graphs of the plasma pressure and magnetic pres- 
sure (top left), the plasma density (top right), the plasma 
temperature (bottom left) and the plasma velocity and 
Alfven speed (bottom right). 

In the trans- Alfvenic case, however, Eq.(|53]) and there- 
fore Eq.(p^) have critical points where M 2 = 1. In order 



to obtain a regular trans- Alfvenic solution it is necessary 
to fix the parameters so that Eq.(53) has no singularity at 
M 2 = 1. The relation 



(A - l)F 2 - 1 



Di 







(55) 



must hold, where a star denotes evaluation of a quan- 
tity at the Alfven point M 3 = 1. Since Eq.(|5f]) involves 
and Pu the parameters Di/G\, D% and A cannot be 
fixed to solve Eq.(|55|) independently of the solution. The 
correct combination of parameter values can be found 
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by a shooting method as seen in Vlahakis & Tsinganos 
(1999). We begin the integration with a reasonable set 
of parameter values and, as the solution approaches the 
Alfven point the solution will asymptote and will approach 
+00 or —00 depending on the sign of the numerator of 
Eq.(|53|)'s right-hand side at this point. Bearing in mind 
that a regular solution will have numerator zero at the 
Alfven point we adjust a chosen parameter iteratively un- 
til we have a regular solution. It is possible to produce a 
regular solution avoiding this iterative procedure by start- 
ing the integration at the Alfven point and imposing fi- 
nite values there, but it is difficult to impose symmetry 
on the loop using this approach. We find it more con- 
venient to begin the integration at the apex, imposing 
symmetry via the initial value of the slope, and iterate to 
find a solution regular at the Alfven point. The parame- 
ter values used in the trans-Alfvenic example given here 
are D x /d = 2.20659184272718, D 2 = 1.0, D 3 = 1.0 and 
A = 2.1. 

Example field lines for the two examples are shown in 
Figs, [l] and ||. The field lines are all loop-like in shape and 
the translational self-similarity of the solutions causes any 
two field lines in a solution to differ from each other only 
by a constant value in the z-direction. The slightly differ- 
ent shape of field line between the two examples is due 
to the different values of D\jC\ and A. In Figs. |3| and |], 
all quantities can be seen to have symmetric profiles. The 
pressure, density and magnetic field strength along the 
field line all have their minima at the apex and increase 
monotonically towards the foot points since these quanti- 
ties all decrease exponentially with height. Note that be- 
cause A > 2 and D2 > the velocity in this example is 
an increasing function of z, as is the Alfven speed (see 
Figs. || and [|, bottom right pictures), and the plasma j3 is 
a decreasing function of z, as well as all being dependent 
on x. This is all in contrast to Tsinganos, Surlantzis & 
Priest (1993)'s example where the velocity, Alfven speed 
and plasma beta are all functions of x alone. (Of course 
the Alfven Mach number depends on x alone both here 
and in Tsinganos, Surlantzis fc Priest (1993).) 

Recall from Subsection 2.2 that non-planar compo- 
nents of the solution are also possible. In the trans- 
Alfvenic example these components are simple: B y must 
be zero while V y is a free function of a and so must be 
constant on a particular field line but can vary arbitrarily 
from field line to field line. In our discussion of trans- 
Alfvenic solutions we confine our attention to the planar 
case By = V y = so that definitions of Alfven points 
M — 1 refer to velocities confined to the x-z plane. In the 
sub-Alfvenic case the absence of critical points permits 
more complicated non-planar magnetic fields and velocity 
fields. The non-planar components B y and V y must satisfy 
Eqs. (||), where 4>a(^4) and L(A) are free functions of a 
which are therefore constant on a particular field line but 
can vary arbitrarily from field line to field line. Examples 
of possible profiles of B y and V y for our sub-Alfvenic ex- 
ample of Figs. and || are given in Fig. ||. In these plots B y 
(left picture) and V y (right picture) are graphed against 



the arclcngth of the largest field line of Fig. [l] and the 
values assigned to the free integrals are $^(^4) = 0.2 and 
L(A) = 0.1. With these values for the free integrals the 
shears of both the magnetic field and the velocity field 
are concentrated near the apex, but B y and V y are not 
directly proportional and so the flow is not field-aligned. 
Note that in any case with or without Alfven points the 
only possibility for field-aligned flow is the planar case 
L = $ A = 0. 

Apart from the existence of critical points in this trans- 
Alfvenic example the examples integrate slightly differ- 
ently from each other. In the sub-Alfvenic case G begins 
from its maximum at the apex of the loop and then tends 
monotonically to 0. In the trans-Alfvenic case, however, 
the non-zero value of D\jC\ and the decreasing value of 
G (and hence of M 2 ) causes the term Di/C\M 2 in the 
numerator of Eq.(|5^) to become large, F' to become pos- 
itive and the loop to loose its loop-like shape. We must 
confine our attention to the region where F' is negative 
if we want a loop-like solution. Adjusting A (and there- 
fore D\jC\) varies the location of this point F' = and 
its value of G(x). This difference in behavior between the 
two examples affects our method of choice of A. 

In this Euler-potential formulation a field line can be 
defined as a line of constant a. We show graphs of quan- 
tities against loop arc length along a single field line for 
each solution and we choose our representative field line as 
follows. In the sub-Alfvenic example where D\jC\ = we 
choose an interval of the x-axis [— Xo,Xo] over which the 
plasma quantities are of a reasonable size, i.e. where M 2 
is not close enough to zero for P or p to be unphysically 
large. In the trans-Alfvenic example where D\jC\ 7^ 
we take the interval [— xo,xo] to be the domain where 
F' is strictly negative. We take as our representative 
field line the line with foot points at the ends of this 
interval x = ±xq. This line is defined by the equation 
a = G(xq) = Gq or equivalently z = hiG(x) — In Go 
(note that the translational self-similarity of the solution 
is explicit in this equation for the field lines). The patterns 
vary slightly from field line to field line but examination of 
several field lines in each example indicates that our cho- 
sen field line is broadly representative in each case, apart 
from the decrease of the plasma (3 with increasing height. 
In the trans-Alfvenic example this variation of (3 limits our 
possible choice of values for A. Because j3 decreases as a 
increases, if our solution is to include a low-/3 field line it 
must extend far enough along the x-axis for small values 
of G{x) (and therefore a) to be included. We find that for 
the present class of solutions A must be taken close to 2.0 
to find field lines with low a and therefore low j3 values. 
This is in sharp contrast to the sub-Alfvenic case where 
F' 7^ away from the apex and a wide range of low-/? field 
lines is available for any value of A. 

Hence in the two examples given here the plasma pres- 
sure is smaller than the magnetic pressure (see Figs. ||and 
^, top left pictures). Given that this is true in particular 
at the apex of the loop where F — and the value of 
G is determined by our wanting a sub- Alfvenic solution 
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Fig. 5. Plots of the y-components of the magnetic field (left picture) and the velocity field (right picture) for the 
sub-Alfvenic example graphed against the arclength of the largest field line of Fig. [[} In this example the free integrals 
are given the values &a(A) = 0.2 and L(A) = 0.1. For both fields the shear is concentrated near the apex although, 
since B y and V y are not directly proportional, the flow is not field-aligned in this (or any non-planar) case. 



or a trans-Alfvenic solution, control (3 here by adjusting 
the value of Pi at the apex as well as adjusting P2 here 
via D2 ■ To keep things simple in this initial study we pre- 
fer to keep the values of the variables at the apex and 
the values of the constants consistent between the two so- 
lutions apart from altering G at the apex to obtain one 
sub-Alfvenic and one trans-Alfvenic example, using the 
low value of A in the trans-Alfvenic and adjusting D\/C\ 
to make the trans-Alfvenic solution regular. Various start 
values of the variables and various values of the constants 
were tried and the results given here are as physically ac- 
ceptable as any found. 

A further difference between the two solutions is that 
their temperature profiles are very dissimilar (see Figs. [| 
and [|, bottom left pictures). Because a is constant along a 
field line and because of Eq. ([so]) the expression for T along 
a field line is of the form KiPi(x)M 2 (x) + K2 where K\ 
and K% are constant. The slightly different Alfven Mach 
number profiles for the two examples (reflected in the 
slightly different field line shapes) and the slightly different 
profiles of the pressure "component" P\{x) have resulted 
in very different temperature functions for the two exam- 
ples. In the sub-Alfvenic example there is a clear tempera- 
ture maximum at the apex of the loop and minima at the 
foot points, while in the trans-Alfvenic example the apex is 
the location of the temperature minimum and the temper- 
ature maxima are close to the foot points. From the inves- 
tigation of many examples it seems that, broadly speaking 
for this special subclass of solutions, temperature maxima 
are characteristic of sub-Alfvenic apexes and temperature 
minima are characteristic of super-Alfvenic apexes. Near 
an apex where M 2 > 1 and F is small, the A term in 
Eq.(^) dominates so that T's slope has the same sign as 
P[ causing T to have a minimum there. However, choosing 



a large value for the constant pressure "component" Pq, 
set to in these examples, affects the temperature profile 
and can even change local minima on the apex to max- 
ima on field lines of small a, although a large value for 
Po affects the plasma (3 especially at large heights. This 
discussion applies only to the special case C2 = of our 
first family of solutions. A more general study of the rela- 
tionship between field line shape and the thermodynamics 
of the more general case C2 7^ and of other solution fam- 
ilies will be the subject of a further paper but these two 
simple examples hint at how much variety is possible. 

4. Conclusions 

We have presented a systematic method for constructing 
exact two-dimensional MHD equilibria with compressible 
flow in Cartesian geometry with a view to modelling struc- 
tures in the solar atmosphere. Using this method we have 
found several new classes of solutions, some of which in- 
clude some well-known solutions as special cases and some 
of which are completely new. As examples illustrating the 
method we have given two simple loop- like solutions, one 
with strictly sub-Alfvenic flow and one with trans-Alfvenic 
flow. The critical points of the trans-Alfvenic example are 
crossed smoothly and all quantities in the model are reg- 
ular. The models have physical properties not previously 
available in exact two-dimensional MHD work. Meanwhile 
two classes of solutions generalise previously-known loop 
and prominence models. It is shown how an analytical 
description of heat input/output may be obtained for a 
solution whether polytropic or non-polytropic. It there- 
fore seems likely that the solution construction method 
may be applied in future in providing quantitative mod- 
els of prominences and arcades to be compared with ob- 
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servations or giving details of the heating profile of the 
structure. Such applications will be the subject of future 
work. 
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